1 COLORS plot

colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines <- c("Balanced" = "solid", "Unbalanced" = "dotted")

2 ESTIMATES VARIANCE REAL DATA

#############################
#PERFORMANCE PARAMETERS VALUE
#############################
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF$Test_environment <- as.factor(data_PERF$Test_environment)
data_PERF$Population <- as.factor(data_PERF$Population)
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] 
data_PERF <- droplevels(data_PERF)
data_PERF$IndicG0<-as.numeric(ifelse (as.character(data_PERF$Gen) == "G0", 1, 0))

nfruit_data <- nlevels(data_PERF$Original_environment)
nhab_data <- nlevels(data_PERF$Test_environment)
npop_per_fruit_data <- nlevels(data_PERF$Population)/nfruit_data
nrep_data <- mean(tapply(data_PERF$Nb_adults,list(data_PERF$Population,
                                                  data_PERF$Test_environment, 
                                                  data_PERF$Generation),length),na.rm = TRUE)
ntrial_data <- mean(data_PERF$Nb_eggs, na.rm = TRUE)
                              

#POISSON VARIANCE VALUE
mperf_poisson <- lme4::lmer(log(Nb_adults+1) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0),
                            data = data_PERF)

sdpop_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[1])                               
sdfruithab_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[3])
sdfruithab_ng_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[2])   
sdfruithab_datapoisson <- sdfruithab_ng_datapoisson
sigma_datapoisson <- sigma(mperf_poisson)


#BINOMIAL VARIANCE VALUE
data_PERF_rate <- data_PERF[data_PERF$Nb_eggs>=data_PERF$Nb_adults,]
mperf_binomial <- lme4::lmer(asin(sqrt(Nb_adults/Nb_eggs)) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0),
                            data = data_PERF_rate)

sdpop_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[1])                               
sdfruithab_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[3])
sdfruithab_ng_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[2])   
sdfruithab_databinomial <- sdfruithab_ng_databinomial
sigma_databinomial <- sigma(mperf_binomial)




###########################
#PREFERENCE PARAMETERS VALUE
###########################
data_PREF <- read.table(file=here::here("data", "DATACOMPLET_PREF.csv"), sep=",", header=TRUE)
data_PREF$Original_environment<-as.factor(data_PREF$Original_environment)
data_PREF <- data_PREF[data_PREF$Original_environment!="WT3",]
dim(data_PREF)
## [1] 4428   12
data_PREF <- data_PREF[data_PREF$Test_environment=="Cranberry"|
                         data_PREF$Test_environment=="Cherry"|
                          data_PREF$Test_environment=="Strawberry",]
data_PREF$Test_environment <- as.factor(data_PREF$Test_environment)
data_PREF$Population <- as.factor(data_PREF$Population)
data_PREF$BoxID <- as.factor(data_PREF$BoxID)
data_PREF$Nb_eggs <- as.numeric(as.character(data_PREF$Nb_eggs))

data_PREF <- droplevels(data_PREF)
data_PREF$IndicG0<-as.numeric(ifelse (as.character(data_PREF$Gen) == "G0", 1, 0))

nfruit_data_box <- nlevels(data_PREF$Original_environment)
nhab_data_box <- nlevels(data_PREF$Test_environment)
npop_per_fruit_data_box <- nlevels(data_PREF$Population)/nfruit_data
nrep_data_box <- mean(tapply(data_PREF$Nb_eggs,list(data_PREF$Population,
                                                  data_PREF$Test_environment, 
                                                  data_PREF$Generation),length),na.rm = TRUE)


#POISSON BOX VARIANCE VALUE
mperf_poisson_box <- lme4::lmer(log(Nb_eggs+1) ~ 1+ (1|Population) + 
                               (1|Original_environment:Test_environment) + 
                               (1|Original_environment:Test_environment:IndicG0) +
                                  (1|BoxID),
                            data = data_PREF)

sdbox_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[1])
sdpop_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[2])                              
sdfruithab_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[4])
sdfruithab_ng_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[3])   
sdfruithab_datapoisson_box <- sdfruithab_ng_datapoisson_box
sigma_datapoisson_box <- sigma(mperf_poisson_box)

#NUMBER OF SIMULS

# Power Analyses
nb_simul <- 5000
#Save the same number of simuls for each SA value
nb_per_SA <- 500
# 
# #Number of simuls for each SA non-genetic value
# tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)

#False positive rate 
nb_simul_fpr <- 100

3 BALANCED AND UNBALANCED DESIGN

3.1 Normal distribution

3.1.1 Balanced

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))

sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))


sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)



## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    3    6   13   16   44   69  100  100  100  100  100  100  100  100  100 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7 
##  100  100  100  100  100  100  100  100  100  100  100   61   26   13    5
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -1.4      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0300000 0.3900000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.5700000 0.5300000 0.6800000 0.6500000 0.6500000 0.6900000 0.6500000 0.7900000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7 
## 0.8000000 0.8600000 0.8300000 0.8688525 0.8846154 0.9230769 0.8000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -1.5      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.1400000 0.1600000 0.2300000 0.2700000 0.3700000 0.3200000 0.5500000 0.5200000 
##         1       1.1       1.2       1.3       1.4       1.5       1.6       1.7 
## 0.5100000 0.6700000 0.6600000 0.6500000 0.7714286 0.5806452 0.8125000 1.0000000 
##       1.8 
## 0.6666667
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -1.5      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.1800000 0.1700000 0.2300000 0.2800000 0.3900000 0.3100000 0.5300000 0.5000000 
##         1       1.1       1.2       1.3       1.4       1.5       1.6       1.7 
## 0.5000000 0.6700000 0.6500000 0.6400000 0.7714286 0.5806452 0.8125000 1.0000000 
##       1.8 
## 0.6666667
tapply(sim$IndicGen, sim$SA_Gen, length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    3    6   13   16   44   69  134  187  315  449  495  520  521  537  560 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7 
##  750 1163 1538 1756 1723 1447 1058  732  477  242  139   61   26   13    5
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
## -1.5 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1 
##    1    1    1    3   15   18   34   75  147  218  292  405  481  530  570  527 
##  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7 
##  533  732 1126 1611 1840 1724 1400 1080  658  486  215  152   70   31   16    5 
##  1.8 
##    3
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -1.4 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2 
##    1    3    6   13   16   44   69  134  187  315  449  495  520  521  537  560 
##  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7 
##  750 1163 1538 1756 1723 1447 1058  732  477  242  139   61   26   13    5
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

3.1.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_normal_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_normal_balancedvsunbalanced

3.1.5 False positive rate analyses

3.1.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho=0, rho_ng=0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.02
mean(sim$IndicGen_aov)
## [1] 0.02
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.01
mean(sim$IndicNonGen_aov)
## [1] 0.02

3.1.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
              rho=-0.1, rho_ng=0,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 
##   3   4  14   6  18  17  13  12   7   4   2
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.2500000 0.0000000 0.5000000 0.2222222 0.3529412 0.4615385 0.3333333 
##         1       1.1       1.2 
## 0.5714286 0.7500000 1.0000000
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.02
mean(sim$IndicNonGen_aov)
## [1] 0.02
tapply(sim$IndicNonGen, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000 
##         1       1.1       1.2 
## 0.0000000 0.2500000 0.0000000
tapply(sim$IndicNonGen_aov, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000 
##         1       1.1       1.2 
## 0.1428571 0.0000000 0.0000000
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##      -0.7      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6666667 0.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##      -0.7      -0.6      -0.5      -0.4      -0.3      -0.2      -0.1         0 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7 
## 0.0000000 0.1111111 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000

3.1.5.3 Non-genetic effects and no genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal",design = "balanced",
              rho=0, rho_ng=-0.1,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.1
mean(sim$IndicGen_aov)
## [1] 0.1
tapply(sim$IndicGen, sim$SA_Gen, mean)
##         -1       -0.9       -0.7       -0.6       -0.5       -0.4       -0.3 
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
##       -0.2       -0.1          0        0.1        0.2        0.3        0.4 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000 
##        0.5        0.6        0.7        0.9        1.2 
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##         -1       -0.9       -0.7       -0.6       -0.5       -0.4       -0.3 
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
##       -0.2       -0.1          0        0.1        0.2        0.3        0.4 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000 
##        0.5        0.6        0.7        0.9        1.2 
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000
tapply(sim$IndicGen, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333 
##        0.9          1        1.1        1.2 
## 0.12500000 0.00000000 0.00000000 0.50000000
tapply(sim$IndicGen_aov, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333 
##        0.9          1        1.1        1.2 
## 0.12500000 0.00000000 0.00000000 0.50000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.33333333 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333 
##        0.9          1        1.1        1.2 
## 0.50000000 0.25000000 0.66666667 1.00000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.00000000 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333 
##        0.9          1        1.1        1.2 
## 0.50000000 0.25000000 0.66666667 1.00000000

3.2 Poisson distribution

3.2.1 Balanced

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))

#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))

#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    1    3    3    6    8    8   15   26   38   46   79   93  100  100  100 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100  100 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
##  100  100  100  100  100  100  100  100  100  100   75   56   41   13   11   11 
##  2.7  2.9    3 
##    6    1    1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -2.2        -2      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0100000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0400000 0.0100000 0.0400000 0.1000000 0.0600000 0.1000000 0.1600000 0.3900000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.4800000 0.5500000 0.7300000 0.7700000 0.8400000 0.9000000 0.9000000 0.9000000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.9200000 0.9400000 0.9733333 0.9821429 1.0000000 1.0000000 1.0000000 1.0000000 
##       2.7       2.9         3 
## 1.0000000 1.0000000 1.0000000
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -2.4      -2.1      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0200000 0.0300000 0.0400000 0.0800000 0.0900000 0.1500000 0.3300000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3700000 0.4600000 0.5700000 0.6300000 0.6300000 0.7300000 0.7800000 0.8400000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.9100000 0.9600000 0.9387755 0.9464286 0.9722222 1.0000000 0.8571429 1.0000000 
##       2.7       2.8       2.9         3 
## 1.0000000 1.0000000 1.0000000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -2.4      -2.1      -1.9      -1.8      -1.7      -1.6      -1.5      -1.4 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -1.3      -1.2      -1.1        -1      -0.9      -0.8      -0.7      -0.6 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0100000 0.0100000 
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0400000 0.0200000 0.0400000 0.0800000 0.1100000 0.1600000 0.3100000 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3900000 0.4400000 0.5400000 0.6500000 0.6400000 0.7300000 0.7900000 0.8400000 
##       1.9         2       2.1       2.2       2.3       2.4       2.5       2.6 
## 0.9200000 0.9500000 0.9387755 0.9464286 0.9722222 0.9629630 0.8571429 1.0000000 
##       2.7       2.8       2.9         3 
## 1.0000000 1.0000000 1.0000000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    1    3    3    6    8    8   15   26   38   46   79   93  117  159  228 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  280  283  316  289  353  303  361  313  326  405  465  645  791  917 1058 1056 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 1079 1062  917  723  609  474  330  287  175  137   75   56   41   13   11   11 
##  2.7  2.9    3 
##    6    1    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
## -2.4 -2.1 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    1    1    1    5    9    9   19   24   33   50   89  100  131  168  201 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  253  277  303  299  349  348  349  320  303  374  454  619  833  985 1080 1088 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 1094  976  869  760  628  444  363  249  162  131   98   56   36   27   14    7 
##  2.7  2.8  2.9    3 
##    4    3    2    1
3.2.1.0.1 TEST
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -2.2   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 -0.7 -0.6 
##    1    1    3    3    6    8    8   15   26   38   46   79   93  117  159  228 
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
##  280  283  316  289  353  303  361  313  326  405  465  645  791  917 1058 1056 
##  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4  2.5  2.6 
## 1079 1062  917  723  609  474  330  287  175  137   75   56   41   13   11   11 
##  2.7  2.9    3 
##    6    1    1
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson

3.2.1.0.2 TEST
################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson

3.2.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim1 <- as.data.frame(t(sim1))

#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)


sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced

3.2.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced_high

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced_high

3.2.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   ggtitle("Poisson data") + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_balancedvsunbalanced

3.2.5 False positive rate analyses

3.2.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho=0, rho_ng=0,   
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.03
mean(sim$IndicGen_aov)
## [1] 0.03
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.01
mean(sim$IndicNonGen_aov)
## [1] 0.01

3.2.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho=-0.1, rho_ng=0,   
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0.04
mean(sim$IndicNonGen_aov)
## [1] 0.04

3.2.5.3 Non-genetic effects and no genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson",design = "balanced",
              rho=0, rho_ng=-0.1,  
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.12
mean(sim$IndicGen_aov)
## [1] 0.12
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4666667 0.4615385 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2000000 0.5000000 0.5000000 0.1666667 0.6666667 0.0000000 0.6666667 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.3       0.4       0.5       0.6       0.7       0.8       0.9         1 
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4000000 0.4615385 
##       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8 
## 0.2000000 0.4000000 0.5000000 0.1666667 0.6666667 0.0000000 1.0000000 1.0000000 
##       1.9         2 
## 1.0000000 1.0000000

3.3 Binomial distribution

3.3.1 Balanced

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))


# Concatenate value
sim <- rbind(sim1,sim2, sim3)


#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   11   79  100  100  100  100  100  100  100  100    7
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0300000 0.4500000 0.6100000 
##       0.3       0.4       0.5       0.6 
## 0.6900000 0.7600000 0.8100000 0.8571429
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
## 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.24 0.33 0.61 0.74 0.75
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##      -0.5      -0.4      -0.3      -0.2      -0.1         0       0.1       0.2 
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0500000 0.3000000 
##       0.3       0.4       0.5       0.6 
## 0.3300000 0.6200000 0.6500000 0.8333333
tapply(sim$IndicGen, sim$SA_Gen, length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   11   79  409 1245 1550 2351 4904 3376  943  124    7
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1    6   75  449 1173 1601 2281 5070 3254  939  139   12
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6 
##    1   11   79  409 1245 1550 2351 4904 3376  943  124    7
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop

#Plot
plot_genetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop

#Plot
plot_nongenetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial

3.3.2 Unbalanced (low)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))

#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 1,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced

3.3.3 Unbalanced (high)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
              design = "unbalanced", disp_design = 4,
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))

# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)




####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA

sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)

#Plot
plot_genetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced_high

################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

#Plot
plot_nongenetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced_high

3.3.4 Plot

data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)


#Plot genetic
plot_genetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_binomial_balancedvsunbalanced

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)

#Plot non-genetic
plot_nongenetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Unbalanced\ndesign") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_binomial_balancedvsunbalanced

3.3.5 False positive rate analyses

3.3.5.1 No genetic and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=0, rho_ng=0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.02
mean(sim$IndicGen_aov)
## [1] 0.02
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0
mean(sim$IndicNonGen_aov)
## [1] 0.02

3.3.5.2 Genetic effects and no non-genetic effects

sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=-0.1, rho_ng=0,
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0500000 0.2926829 0.4062500 0.5714286
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0500000 0.2926829 0.4062500 0.5714286
## Local adaptation (non-genetic)
mean(sim$IndicNonGen)
## [1] 0
mean(sim$IndicNonGen_aov)
## [1] 0

3.3.5.3 non-genetic effects and no genetic effects

### 
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
              rho=0, rho_ng=-0.1,
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_databinomial,
              sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial, 
              sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.07
mean(sim$IndicGen_aov)
## [1] 0.07
## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0000000 0.2075472 0.3333333 0.6000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.1       0.2       0.3       0.4 
## 0.0000000 0.3018868 0.4074074 0.6000000

4 PLOT Balanced vs Unbalanced

POWER_ALL<-cowplot::plot_grid(plot_genetic_normal_balancedvsunbalanced,
                              plot_genetic_poisson_balancedvsunbalanced,
                              plot_genetic_binomial_balancedvsunbalanced,
                              plot_nongenetic_normal_balancedvsunbalanced,
                              plot_nongenetic_poisson_balancedvsunbalanced,
                              plot_nongenetic_binomial_balancedvsunbalanced,
                              ncol = 3, nrow = 2,
                              scale = c(1,  1))
POWER_ALL

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_BALANCED_and_UNBALANCED",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   POWER_ALL, base_height = 20/cm(1),
#                   base_width = 40/cm(1), dpi = 1200)

5 BOX EFFECT Absent vs Present

5.1 Non box

rm(sim)
rm(data_proportion_gen, data_proportion_nongen, 
   sim_select_genetic, sim_select_nongenetic,
   data_prop, data_prop_gen, data_prop_nongen)

#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, 
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))


# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines_box <- c("Present" = "dotted", "Absent" = "solid")
####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0
data_prop_gen <- data_prop



################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0
data_prop_nongen <- data_prop

5.2 Box (low)

#With lower value of sdbox (close to sdpop): 
sdbox_datapoisson_box
##     BoxID 
## 0.4464253
sigma_datapoisson
## [1] 0.8727978
sdpop_datapoisson
## Population 
##   0.504277
sdfruithab_datapoisson
## Original_environment:Test_environment:IndicG0 
##                                     0.8111606
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = sdbox_datapoisson_box,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)


####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0.5
data_prop_gen <- rbind(data_prop_gen,data_prop)


################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

5.3 Box (high)

#With lower value of sdbox (close to sdpop): 
sdbox_datapoisson_box
##     BoxID 
## 0.4464253
sigma_datapoisson
## [1] 0.8727978
sdpop_datapoisson
## Population 
##   0.504277
sdfruithab_datapoisson
## Original_environment:Test_environment:IndicG0 
##                                     0.8111606
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.5, rho_ng = -0.5, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = 0, rho_ng = 0, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data, 
              nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
              sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson, 
              sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)

sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)


####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))

data_prop$sd_box <- 1.5
data_prop_gen <- rbind(data_prop_gen,data_prop)



################################# 
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 1.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)

5.4 Plot

#################################
### Common plot
#################################
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$sd_box <- as.factor(data_proportion_gen$sd_box)


#Plot genetic
plot_genetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_gen, 
   aes(x = val_SA_genet, y = proportion, color = method, linetype = sd_box)) + 
   geom_point(size = 2) + 
   geom_line(size = 1, alpha = 0.7) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Box effect\n(sd box)") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_ALL_box

#Data non-genetic 
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method, 
                                      c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$sd_box <- as.factor(data_proportion_nongen$sd_box)



#Plot non-genetic
plot_nongenetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_nongen, 
   aes(x = val_SA_nongenet, y = proportion, color = method, linetype = sd_box)) + 
   geom_point(size = 2) + 
   geom_line(size = 1) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Box effect\n(sd box)") +
   scale_color_manual(values = c("#459B76","#DCA237")) +
   scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_ALL_box

power_poisson_box<-cowplot::plot_grid(plot_genetic_poisson_ALL_box,
                                      plot_nongenetic_poisson_ALL_box,
                          labels=c("A", "B"), 
                          label_size = 15,
                          ncol =1, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_poisson_box

#Save plot
name_plot<-paste0("Simul_powertest_poisson_boxeffect_",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_poisson_box, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)

6 SIMUL WITH REAL UNBALANCED DESIGN

6.1 Dataset

#unbalanced dataset
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] #mvrnom doesn't work with unbalanced matrix
data_PERF <- droplevels(data_PERF)
tapply(data_PERF$Population, list(data_PERF$Generation, data_PERF$Population), length)
##    Blackberry31 Blackberry32 Blackberry33 Blackberry34 Blackberry35
## G0           38           42           25            4           NA
## G2           27           39           44           11           46
##    Blackberry36 Blackberry37 Blackberry38 Blackberry39 Blackberry40
## G0            8           15           21           10           15
## G2           17           33            5           24           30
##    Blackberry43 Blackberry44 Blackberry45 Cherry103 Cherry104 Cherry3 Cherry47
## G0            8           12           43         3        39       9       36
## G2           10           30            4        19        16       9       39
##    Cherry50 Cherry51 Cherry52 Cherry6 Cherry7 Strawberry42 Strawberry44
## G0       36       48       45      36       6           15           21
## G2       11        1        4       7      37           46           30
##    Strawberry53
## G0            6
## G2           63

6.2 Normal distribution

####################################
###### Simulation
####################################
#Simul data 
sim <-  sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "normal", 
               unbalanced_dataset = data_PERF, 
               sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, 
              rho = -0.05, rho_ng = -0.05, 
              sigma = 0.5)
             

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 
##   2  58 240 568 819 944 827 675 421 240 127  50  21   8
#Save the same number of simuls for each SA value
nb_per_SA = 50
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                     aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data (unbalanced experimental design)") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 
##   1  60 236 528 890 969 826 617 430 215 143  47  26   8   4
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


#Plot
plot_nongenetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_normal_unbalanced_real

6.3 Poisson distribution

####################################
###### Simulation
####################################
#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "poisson", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 1, sdfruithab_ng = 1, rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 2.2 
##   3  26  66 123 195 289 353 422 461 466 451 394 387 347 253 233 166 109  86  71 
## 2.3 2.4 2.5 2.6 2.7 2.8 2.9 
##  34  29  15   9   6   4   2
#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                      aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_genetic_poisson_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 
##   1   4  29  61 123 174 256 386 449 495 478 443 405 379 326 239 201 174 107  93 
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9   3 
##  63  40  28  17  14   7   2   3   3
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_poisson_unbalanced_real

6.4 Binomial distribution

####################################
###### Simulation
####################################
#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "binomial", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, 
              rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 
##   2  58 240 568 819 944 827 675 421 240 127  50  21   8
#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                  aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors) 
plot_genetic_Binomial_unbalanced_real

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 
##   1  60 236 528 890 969 826 617 430 215 143  47  26   8   4
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 50
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_Binomial_unbalanced_real

6.5 All plots unbalanced

power_plot_all_unbalanced_real<-cowplot::plot_grid(plot_nongenetic_normal_unbalanced_real,
                                   plot_nongenetic_poisson_unbalanced_real,
                                   plot_nongenetic_Binomial_unbalanced_real,
                                   plot_genetic_normal_unbalanced_real,
                                   plot_genetic_poisson_unbalanced_real,
                                   plot_genetic_Binomial_unbalanced_real,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_all_unbalanced_real

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_UNBALANCED_real",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_all_unbalanced_real, base_height = 20/cm(1),
#                   base_width = 42/cm(1), dpi = 1200)
#